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Abstract. Typical features of the Transmission Line Matrix (TLM) algorithm 
in connection with stub loading techniques and prone to be hidden by common 
frequency domain formulations are elucidated within the propagator approach. 
In particular, the latter reflects properly the perturbative character of the TLM 
scheme and its relation to gauge field models. Internal 'gauge' degrees of freedom 
are made explicit in the frequency domain by introducing the complex nodal S- 
matrix as a function of operators that act on external or internal fields or virtually 
couple the two. As a main benefit, many techniques and results gained in the time 
domain thus generalize straight away. The recently developed deflection method 
for algorithm synthesis, which is extended in this paper, or the non-orthogonal 
node approximating Maxwell's equations, for instance, become so at once available 
in the frequency domain. In view of applications in computational plasma physics, 
the TLM model of a relativistic charged particle current coupled to the Maxwell 
field is traced out as a prototype. MSC-classes: 65C20, 65M06, 76D05 
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1. Introduction 

When P. B. Johns and coworkers introduced the transmission line matrix 
(TLM) numerical method in the early seventies [7] they were certainly not 
aware of handling gauge field concepts. In fact, the perturbative aggrega- 
tion of internal gauge degrees of freedom plays a salient role in TLM in the 
form of 'stub loading'. This well known technique, ad hoc applied in a great 
many variants to model continuous mesh deformations or material param- 
eters, e.g., expands to a powerful framework that allows for a systematic 
integration of manifold interactions. Stubs in a generalized sense may for 
instance add mass to a field much alike a Higgs field [15] , or induce super- 
conducting or non-reciprocal behaviour [13) , [14j . Very much of the strength 
and versatility of TLM rests on this feature that indeed relates the method 
to some pioneering developments in modern particle physics. 

Varied formal characterizations of the TLM scheme have been given since 
the early transmission line network interpretations. Some authors proceed 
from Huygens' principle 8J (sometimes misunderstood) others use concepts 
of state space control theory [TO] or random walk [2] , [PJ • The present paper 
is not intended to add a further characterization to these and alternate 
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interpretations which are usually built on the cinematics of a special kind 
of propagating quantities (waves, particles, temperatures etc.). Instead, this 
study focusses on the internal structure of the TLM algorithm irrespective 
of its physical interpretation. It nevertheless remains remarkable that so 
many heterogeneous perceptions exist and apparently match the algorithmic 
scheme of the method. The perturbative structure of TLM, with its gauge 
fixture not at least, certainly explains much of that multiplicity. 

With sustained evolution of computer performance in coming years there 
will be an increasing need of highly flexible numerical schemes, complemen- 
tary to powerful standard solvers. The versatile structure and the diver- 
sity of yet known applications in wave propagation [1] , chemical reaction- 
diffusion [TO], heat transfer [5], plasma physics [T5], and further domains 
provide some evidence that the TLM method can play a competent role on 
the court. 

This requires, however, a refined technical outfit being built on sound 
mathematics. Theorems like the deflection lemma [14], a generalized ver- 
sion of which is given in this paper, should constitute a canonical basis for 
algorithm synthesis. The usage and usefulness of those tools amenable to 
the propagator approach have already been shown up within a series 
of numerical applications [T^] - [T5]. On the issue of this paper the TLM 
model of a relativistic charged particle current coupled to the Maxwell field 
provides supplementary illustration. 

2. The algorithm and its interpretation 

This section is widely heuristic and starts slightly outside our subject. Con- 
ceptually, the TLM algorithm will be separated from its interpretation in 
terms of any modelled fields. This obviously implies outlining a canonical 
passage between two severed parts: physical fields and computed quantities. 
The TLM typical manner of generating these quantities, i.e. the computa- 
tional part, is then headed for necessarily on a level of abstraction beyond 
the commonly considered Maxwell field. 

It is in the obvious nature of any splittings (conceptual just as physical) 
that junctions appear double faced, and so the joints between the TLM 
algorithm and its interpretation, viz. ports and nodes in cells, take a twofold 
role. So it has to be clarified what is meant in either context. 

At the interface to physical fields in space a cell stands for a convex 
polyhedron in configuration space M 13 and a port for a distribution, i.e. a 
continuous linear functional acting on a field, with its support on a closed 
subset of the cell boundary (a point, line segment, or face e.g.). The dis- 
tributional support is also referred to as the location of that port. If no 
confusion can arise it will be simply identified with the port. 

Again, in the same context, the node comprises the family of all port dis- 
tributions of a cell, yet shifted by translations into its interior (where concep- 
tually one effective internal state of the cell at a time interacts with bound- 
ary states - this touches the algorithmic role of ports and nodes treated 
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below). Note that in many cases distributed integrals instead of point sup- 
ported functionals, viz. the shifted nodal images of the ports, probe the 
fields in the cells. 

In a condensed slightly formalized writing: Every port p of a cell together 
with the node gives rise to a pair of linear forms (p, Z), (p, Z)~ which probe a 
real or complex vector field Z = (Z\, . . . , Z m ); m € N+, at the cell boundary 
and within the cell. The two forms are naturally interrelated by pull-back 

(i) ( P ,zr = ( P o S ,z) = (p.Zor 1 ) 

where s denotes the pertinent translational shift of the port onto its nodal 
image, cf. fig 1. Defining these forms consistently as distributed or point 
supported functionals (finite integrals, e.g.) is a constituent part of any 
interpretation of the TLM model. 



node 



cell face 




Figure 1: Cell boundary ports and their nodal images 

Quite different, and still more abstract, are the roles of ports, nodes, 
and cells in the computational context. Every cell stands here for a direct 
sum (viz. a triple) of linear spaces, which respectively represent observable 
link states in the cell or on its boundary, and unobservable gauge (or stub) 
degrees of freedom. Observable vectors z v ^ n are labelled by ports p or their 
nodal images n, and are linearly interrelated to physical fields by an inter- 
pretation. In a more precise diction, an interpretation may be defined as 
a family of linear maps Iq, indexed by the cells of a mesh, with domain a 
distinguished class of vector fields in configuration space which represent 
physical fields and with their range in the set of observable cell states: 

(2) h : Z, - ; z£ = ( p„ Z,) , 2 « = ( p„ 

(jU labels the ports and pertinent fields in cell £ and superscripts p, n refer 
respectively to ports on the cell boundary or to their nodal images in the 
cell - a port being counted here as many times as it acts on a different field.) 
Note that no fixed linear relations of such a kind interconnect any gauge 
vectors with fields. 

Ultimately, and by it most markedly characterizing the TLM scheme, 
the ports together with their nodal images represent scattering channels 
for states incident at the cell boundary and scattered in the node. As a 
scattering channel every common port of two adjacent cells thus naturally 
connects their nodes, i.e. constitutes a link. 

The inherent scattering scheme of the algorithm postulates in the line 
of Lax-Phillips [3] a direct splitting of the link states into incoming and 
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outgoing components 



(3) 




which of course has to be defined in harmony with the propagation proper- 
ties of the modelled fields. In reference [14] the details are carried through 
for the Maxwell field, for example. The Poynting vector evaluated at the cell 
faces yields there a canonical splitting and a quadratic form that is positive 
or negative, respectively, on the incoming or outgoing subspaces. 

In TLM diction the latter are usually called the (spaces of) incident 
and reflected quantities and in the following 7Tj ra , 7r out denote the pertinent 
projections. 

The TLM scattering algorithm, which is eo ipso perturbative, solves 
finite difference equations in the form of model equations between states 
that in any interpretation should approximate physical fields (i.e. solutions 
of integro-differential equations) in any sense and order (sect. 4). This is the 
requirement of consistence, as a cardinal property of the interpretation. 

In contrast, the question of internal convergence, which is essentially 
that of computational stability, refers to intrinsic properties of the TLM 
algorithm such as are subject of this paper. That question has thus to be 
taken up in due course. (We are not unhappy that it is in general easier.) 

3. The propagator relations 

Causality requires in the time domain the scattered states of a cell to be 
functions of previously incident states. With a cell dependent reflection op- 
erator 31 the outgoing port quantities are thus 



where r denotes the time step. (Suitable definitions release from introducing 
line impedances, such as are commonly used instead of the projections ([3]) 
- in fact they are substitutes for projections, cf. the appendix.) 51 assembles 
the physical interactions within the cell in a geometry dependent manner, 
and represents in this sense the TLM system locally, while the global topo- 
logical structure, as given by the port linking scheme and the boundary 
conditions, are brought in by the connection map C. At every time step 
the latter transmits the reflected port quantities to adjacent cells or, at a 
system boundary, back into the same cell and contingently adds the field 
excitations by acting linearly as 



where £, r\, i? label the cells and z\ xc denotes any excitation. 

Short-cut, every pair (C,3?^) consisting of a connection map C and a 
family 31^ of reflection operators labelled by the cells of the same mesh will 
be referred to as a TLM system. 



(4) 



*Lt (t + T)=H(zf n (t-HT)) 



(5) 
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The interesting objects in a TLM system are, of course, the causal op- 
erators which like the pages of a book are merely bound together by 
the 'book binding' function C. In addition to being defined on back-in-time 
running sequences of incident states, 3? may itself depend on time and may 
be non-linear. Before venturing on the TLM typical structure of 3? a glance 
is cast at some cinematical relations between port and node quantities. 

Operationally, the TLM approach, as a genuine scattering scheme, splits 
off a cinematical part from the dynamical evolution by handling some inter- 
actions as a perturbation. In classical Lax-Phillips scattering, for instance, 
two 'intertwingled' one parameter unitary groups describe respectively the 
free and interacting dynamics [3]. A related ansatz leads in perturbative 
quantum field theory to the Dyson expansion of S-parameters in the form 
of Feynman diagrams which again couple free fields [4]. Technically, the 
TLM algorithm operates similarly with freely propagating port and node 
quantities that before and past reflection are forced into fixed phase rela- 
tions by imposing in the time domain 

(6) > + 

Note well that these are technical arrangements to enable the scattering 
scheme, while only total quantities, 

(7) *"' n = *£ n e*£3 , cf. © 

represent observable states and in fact enter the model equations (cf. sect. 4). 

Scarce and innocent as they appear, settings ((4J El [6j) together with ((3J 
contain yet very much of the substance of TLM. On the basis of essentially 
these ingredients, which form the core issues of the so-called propagator ap- 
proach [11], [15], internal 'gauge' degrees of freedom arise naturally on intro- 
ducing additional interactions by a perturbation, along with S-parameters 
that couple these to the original states. A simple way of giving a rigorous 
meaning to that is the deflection lemma. Originally introduced in reference 
P3], the lemma is subsequently generalized to handle higher order equa- 
tions. 

4. The deflection lemma 

Finite difference equations that in the form of model equations approximate 
the dynamical equations of any physical interpretation interrelate in the 
time domain sequences of total node and port quantities 

(8) [^• P ] = (^ ,P (i-Mr)) M=0ll , 2) ... 

which represent observable states in the cell and on its boundary, at present 
time t and in their history. A maximum /i that eventually enters the equa- 
tions obviously reflects the temporal order of the modelled dynamical prob- 
lem, which may be finite or infinite. To any order the model equations can 
be written 



(9) 



3K][^] = o 
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with a suitable causal operator 5" (where causal stands always for being 
defined on back- in-time running sequences such as l[8])), subscript +(— ) 
denotes a r/2 time shift 

(10) K f_ ,] := [^] t+( _ )f = (t+(-)I - M r))^ 

which in harmony with the phase relations © synchronizes port and node 
switching in equations ©, for every vector zf n (t) incident at a cell and 
constant on each time interval [kr, (k + l)r); fc S Z. Note that 5" may itself 
be time dependent in like manner, i.e. as a time step function matching 
integer multiples of r. 

In the literature equations ((9|) are encountered in manyfold guises. For 
instance, in reference [IT] Maxwell's integral equations discretized in a con- 
vex hexahedral cell of else arbitrary shape take the form 

(11) ipozP(t)=<poz n (t + ^) +viz n (t-^) 

where z — (u, i) represent finite integrals over electric and magnetic field 
strengths, "00 essentially interchanges u and i and ifo, ipo are selfadjoint 
M-linear operators, cf. references [S],[TS] and the appendix of this paper. 

By saying that a TLM system (C, Dlf) generates solutions of equations 
([9]) it is meant that the total fields z n ' p , cf. (|3|7|) . satisfy the equations for 
all vectors zf n incident at any cell (and as step functions of time of course 
matching r) . Obviously, this defines a local, so far purely algebraic property 
of the TLM system, i.e. one referring to 3^ and indepedent of C that still 
disregards all questions of convergence or computational stability, for in- 
stance. As a rule, the TLM system solving ((9J depends strongly on the time 
step and additional contraction properties must be checked on selecting a 
particular stable system (cf. sect. 5). 

The deflection lemma establishes general recurrence relations between 
solutions of model equations that differ by a perturbation. 

Let a TLM system (C,3^) generate solutions of equation ([9]) with any 
9" of the form 

(12) [z*>] = ]T <p»z n ( t + 1 - M r ) + (t - M r) 

MSN 

wherein (p„, ip^ are K-linear operators (many if not almost all of which may 
be zero) . Consider then a perturbation of equation © by a causal possibly 
non-linear operator 3 that maps into the image space of 5" (i.e. y>„, ip^ and 
3 share the same image space). 3 and 3 may be time dependent, i.e. step 
functions of time matching r. 

Then one can state the following 

Proposition (Deflection Lemma) A TLM system (6,3^), with the 
same connection map as above, generates solutions of equation 



(13) 



3\z»][z p ]=3[z1][z p ] 
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if and only if the deflection 



(14) 



D := X 



satisfies recursively 



(15) 



t — flT 



Corollary. Let </?o be injective (one to one) on the outgoing subspace and 
assume ipo o ir out covers the ranges of 3, ° ^out, and ip^ o Tr out ; fi G N. 
Then 



(16) ■D t -=-^ o7T out ); 1 {2t[z!][z p ] + &n+^-i) t Vt-^} 



is the reflection map of a TLM system (C, 31^) which generates solutions of 
equation (fT5|) . 

Moreover, if Do = 0, then this system is the unique one which branches 
from (6,3^) at t = 0, i.e. that initially coincedes with the unperturbed 
system. 

Condensed statements of that kind support a comment before being 
proved: 

Guided by conceivable applications the perturbation 3 is written as a 
function of total port and node states and of time. Most commenly, it is a 
function of nodal states only. 

If ipo o Tt out is invertible, then the negative time shift in <| 1 3|) makes 3t 
depend on nodal states up to 



thus ensuring that D enters the right-hand side of recursion (|16p (implicitely 
in the argument of 3 t ) up to CD t _ r a t most, rather than up to D t which ob- 
viously would lead to ill-defined recurrence relations without further severe 
restrictions on 3- 

Note that from K-linearity only additivity of tp^ , ip^ enters the following 

Proof Let hence (C,$£) generate solutions of (J9j) and let for a second TLM 
system (C, 3?^) the deflection D be defined by (Q3|. Then for every sequence 
of states [zf n ] incident at a cell the total states of the second system are 



defines recursively a causal operator D such that 



:= ft + D 



(18) 




t-T 
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This substituted for z n ' p in equation (fT3|) yields 

( 19 ) = E M£N { M* n (* + i-/*r) + l»t- Mr ) 

+V M (2 p (t - /*t) + D t _ T _^-)} 
= 3t[zl][z"] 
By virtue of the additivity of tp^, \n follows 

(20) 2[z»] [F] = S\zl] [z p ] + J2 f^t-^r + ^Vt-T-vr = 3t [zl] [z p ] 

Equation © states that ^[z™] [z p ] above vanishes and inspection of the 
remaining terms yields the proposition. 

5. The S-matrix propagator 

The first order linear equation (jlip admits a solution in closed form as an 
iterated scattering process. This has been pointed out in detail in references 
03] > [S] a t the example of the non-orthogonal Maxwell field model. Using 
the commuting family of projections which split the state space of each cell 
Q, with identity operator Id^, into observable link states and unobservable 
gauge or stub vectors 

(21) id c = n © tt s 

and then again the link states into incoming and outgoing subspaces of port 
and node quantities 

(22) 7T; = TT in © TT out 

the reflection map of the first order linear process solving is 

z P ou& + t) =*[*?„] = z n out {t+ T -) 
= n out S z" n (t+ - ) 

OO 

(23) +7rout sY,(KsSTzZ l (t+~-HT) 

S denotes the nodal scattering operator which acts between in and outgoing 
node and stub vectors and can thus be written 



(24) S = {n out © tt s ) S (n ln © tt s ) = K + L + M + N 



Time and frequency domain TLM 



where as usual 

(25) K := n out SiTi n , L := ir out SiT s , M := ir s STr in , N := n s Sn s 

have been introduced. With these operators equation (I23p reads more con- 
cisely 

oo 

(26) Ct (*) - Kz^ (t)+Lj2 N"Mz? n (t - ( M + 1) r) 

fj,=0 

(graphically represented in fig. 2a). 

K, L, M, N are easily determined from equation (fTT]) by substituting 
for z n ' p = Zfo P + z out m this equation the total scattering response of a 
Dirac pulse 

t \ Zn i f t £ \o. t) 

(27) i(t)=C(i + ^) = [ ' > 

1 [ else 

applied and added to (|2"51 |2"1)|) . By evaluating equations (JTTJ) with these 
replacements at time intervals [fer, (fe + 1)t), fe s N, provides operator iden- 
tities that are independent for k = 0, 1 and k > 2 under rather general 
conditions, cf. [12j - [15] and which allow to determine K uniquely, while 
L,M, and iV are unique at most up to a linear automorphism of the stub 
space. Indeed, sheer inspection of the scattering response (f2"51 12"6")) makes 
clear that any invertible linear transformation G of the stub space yields an 
equivalent scattering representation of H if S is replaced by 

(28) S' := K + L' + W + N' 

with V := LG~ X , M' := GM, N' := GNG^ 1 . It is this property that led 
to the understanding of stub vectors as internal gauge degrees of freedom, 
the automorphisms of the stub linear space (in transmission line parlance: a 
group of generalized impedance transformations) constituting the pertinent 
gauge group. 

As already mentioned, the inherent gauge features of the TLM scheme 
can actually be traced back to the deflection lemma, the latter displaying in 
like manner the generation mechanism for stub degrees of freedom as rather 
suggestively also their gauge property. 

Given a TLM system (6, JIq) that generates solutions of equations ((9]) 
with any operator tpa sharing the technical prerequisites refered to, it is 
a straightforward exercise applying the lemma and its corollary to write 
down the updating instructions for a 'deflected' system (C, that solves 
the perturbed equations (fl3|) : 

zl ut (t + r) :=% [<]+D(t) 

(29) I(t + r) := - fa o ir out );^ T 3 t+T [^][z p ] 

D (t + t) := / (t + t) - (<p o Tr out );+ T ^ + ¥V+i) t+r D (* ~ I*t) 
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with all updating operations done in this order and implicitly after the first 
line 

Z n = 4a (*) + ^out (t + T) 

Besides the proper interaction term /, which essentially appears in the model 
equations () 1 3(1 and is thus amenable to a direct physical interpretation, the 
deflecting field D introduces additional degrees of freedom that in general 
remain without counterparts in the equations. 

The gauge character of the deflection field is suggestively illustrated 
within scattering representations (|23l26p for H and - assuming for once 
that such exist. If 5~ denotes the S-matrix of then it admits a block 
decomposition in matrix form 



(30) 



s 









wherein S represents the S-matrix of 31 and N~ act on the additional 
deflection field D. The updating instructions ([29]) obviously fit with L~ of 
the form 



(31) 



1 0. 



•°7 



(all other matrix elements of L~ being zero). Any invertible transformation 
G, according to (f28|) simultaneously applied to L~, M~ and N~, then 
cannot alter and thus yields an equivalent representation. (The reader 
may trace this forth to a set of modified updating equations (|29p .) 

Spatial deformations of the TLM mesh continously get through stubs. 
In fact they are modelled using stubs [T] , [5] , [T^] . Therefore smooth connec- 
tions between mesh point dislocations and infinitesimal impedance trans- 
formations G exist that textually complete the gauge field analogy in the 
sense of fibrations [5] . 

Once the gauge fixture of TLM being made available, full advantage can 
be taken of its power and flexibility. Aside from the nonlinear potentiality 
opened by the deflection lemma, manifold variants of the linear scheme (|26l) 
are conceivable. For instance the generalized ansatz 

fc 

(32) z n out (t) = KzZ (t) + E L - E N £ M *tfn (* - 0" + K ) T ) 

k=i /ieN 
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(which is graphically depicted in fig. 2b) is suited to solve linear model equa- 
tions of the class (fl2|) with time differences in the port states up to order k. 



N 
T 

M * L 3x 2 V 

*■ v 2 < 

• - Kut «?„ - • r 

T 
K 



a) b) 
Figure 2: Linear scattering schemes 

Every path from left to right following the arrows yields a delayed or phase 
shifted contribution to the scattering response. (Loops may be repeatedly run 
through any number of times.) 

A direct way to prove the above statement and to derive S-matrix blocks 
K, L K , M K and N K in (f3"2")l is again by evaluating equation (TT2"]) at successive 
time steps for an incident Dirac pulse (|27p . 

6. The stable process in the time and frequency domain 

Parallel to the TLM typical splitting of states into link and gauge degrees of 
freedom external and internal conditions for computational stability may be 
distinguished, which respectively stress some physical aspects and a purely 
numerical origin of instability. 

In any interpretation the total link vectors are directly related to phys- 
ical quantities that are in general subjugate to conservations laws. It is 
hence merely a strengthened consistence postulate to require that the re- 
flection map of each cell must comply with these laws. The entire reflection 
response of an incident Dirac pulse integrated over time has for instance to 
be physically balanced with the input. 

Actually this goes along with the Lax-Phillips contraction conditions 
[3] which are sufficient to prevent the observable fields from increasing in- 
definitely in the mesh. To implement these conditions numerically in the 
TLM spccifical fashion, internal algorithm stability has in addition to pre- 
vent gauge fields (i.e. the 'stub quantities') from piling up to infinity. For 
instance, in a stable non-linear perturbed process which solves equation () 13(1 
the deflection must die out in time for any transitory incident fields that do 
not pass over a certain level. 

A universal condition that guarantees computational stability, even in 
the non- linear case (at least up to some maximum level of incidence), is quite 
evidently all updating operations must be contractive in a neighbourhood of 
zero. For the linear processes (|2"tj)) . (j3"2")l the operators N( K ) control updating 
between gauge fields. Hence a sufficient condition for internal stability of 
these processes is 

(33) || N ||< 1 
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if the norm |[ . . . || is submutiplicative, i.e. characterized by 

(34) || N 1 - 1 || < || N f , n e N 
The sup-norm with respect to any state space norm || z || 

(35) ||JV|| BUp := sup || JVaf || 

11^= ll=i 

or the Hilbert (spectral) norm 

|| N \\h '■= max | A 

A eigenvalue of N 



are, for instance, submultiplicative in the sense of (|34|) . All operator norms 
occurring in the following are tacitely assumed to share this property. 

Condition (|33|) is in general realized in the cases of interest by intro- 
ducing bounds for the time step. In the frequency domain, for a process of 
angular frequency ui = 2tt/, i.e. 

(36) z^ out (t + ht) = e^z^ out (t) , a E Z 

equation (|26|) reads 



(37) z n out =(K + LJ2 N»Me^+^ z » 

V p=0 / 

With C-linear operators K, L, M, N this turns into 

(38) Ct =U+ e-^L £ (e-^Nf M j 

A stable process, characterized by condition (|33p . ensures that the power 
series in e~i u ' r ~N converges to (Id—e^^^N)^ 1 . Equations (f3"5| thus becomes 
Zout = S~z? n with 

(39) S~ =K + L (e juJT Id - AT) _1 M 

Evidently, in the frequency domain the complex S-matrix connects directly 
incident with outgoing fields and stubs do no longer enter the algorithm 
in the form of explicitely computed quantities. Instead, the stubs yield im- 
plicitely an additive contribution to the S-matrix 

(40) S~ = L (e JUT Id -Ny 1 M 

extracted from decomposition ()39[1 . It is virtually this gauge contribution 
that anew characterizes the TLM scheme here in the frequency domain. 

Of course, knowing the nodal S-parameters of each cell, 'intrinsic S- 
matrices' of entire mesh systems can be computed, just as then S-parameters 
for any combination of submeshes. (Proceeding in that direction opens 
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rather a special field of application to linear network theory than deeper 
structural insight into the TLM method). 

In the frequency domain the TLM algorithm solves the eigenvector equa- 
tion 

(41) (es c ~v WT - id) ( Zin ) c = o 

where C runs over all mesh cells and the connection map C interlaces a steady 
state excitation which may be imposed on the mesh boundary according to 
any mode templates for instance. 

Various methods for solving equations (|41|) are imaginable. We experi- 
mented with classical conjugate gradient schemes, such as Fletcher-Reeves 
and Polak-Ribiere which, however, showed poor convergence for larger mesh 
systems of several thousand cells. Very satisfactory results have been ob- 
tained, in contrast, with a fixed point iteration analogous to a static time 
domain process (of course with complex fields). Indeed, reiterating the in- 
structions 

{Zin)^ ^ C (z ou f -\- Zexc) ^fc 

(42) (*a»t) c ,*+i = e jW S c ~ (z m ) Ck 

converges for any steady state excitation (z exc )^ quite rapidly to a fixed 
point, and thus to a solution of equation (|4"T|) which matches the imposed 
excitation pattern. Still faster convergence had been attained by starting 
with zero fields and smoothly approaching the steady state excitation over 
a transitory phase of some hundred iterations. 

Fig. 3 displays convergence thus attained for a 3 dB waveguide directional 
Riblet coupler. 



Riblet coupler: sll at steady state 14.25 GHz [measured: -41.6 dB] 
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Riblet coupler: s21 at steady state 14.25 GHz [measured: -3.08 dB] 
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Riblet coupler: s31 at steady state 14.25 GHz [measured: -3.12 dB] 
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Riblet coupler: s41 at steady state 14.25 GHz [measured: -38.4 dB] 
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Figure 3: Convergence in the frequency domain of a fixed point iteration |42J 
for a 3dB waveguide directional coupler (Riblet coupler). The deviations 
from measurements at the real coupler are mainly due losses, which are 
neglected in the computational model. 

7. Charged particles coupled to the Maxwell field 

It has been stressed that this study focusses on the intrinsic structure of the 
TLM algorithm, not on any interpretation. Still one application is sketched 
in this section to illustrate the excellent capacity of TLM to cooperate with 
a tied finite difference scheme. Beyond that, a variety of original research 
applies concepts of this paper and should be consulted for technical illustra- 
tion. In particular, part of references [llj - [15] constitute by now 'standard 
technicality' which can be referred to in the following. (A bridge between 
the non-orthogonal Maxwell field model, as treated in [15] , and the presently 
generalized formalism is built in the appendix.) 

The following application extends the Maxwell field model in a non- 
linear fashion. By deflection, a relativistic charged particle current is coupled 
to the apriori free electromagnetic field. The underlying dynamical relations 
are thus Maxwell's equations and the equations of motion for a particle 
current with interactions being brought in by Ampere's law and the Lorentz 
force. 

The relativistic mass of a particle at velocity v is 



(43) 



m ( 1 



)' 
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where mo and cq respectively denote the particle rest mass and the velocity 
of light. For particles of charge qo the Lorentz force in the laboratory frame 
is 

(44) F L = j t (mv) = q Q (E + v A B) - v c mv 

where a mean collision frequency v c has been introduced to take internal 
friction into account; cf. [Sj. 

Equation (|44[) is linear in dv/dt and reads explicitely 

d , dv dm „ , dv „ 

(45) * (mv) = m ^ +V ^ M ^= Fi 
with 

(46) M = m(6ij + \v i Vj) i j =lt2i3 ; A = (c{) - v 2 ) 1 

by virtue of (|43p. M is obviously selfadjoint and strictly positive for v < cq, 
hence invertible, and yields the relativistic particle acceleration as 

(47) v = M -1 Fz 

Then n particles per unit volume of (mean) velocity v generate a current 
density 

(48) j c = p ■ v ; p = nq 

which has to be included in the generalized Ampere's law (1st Maxwell's 
integral equation) 

(49) J lids = J £^dS + J K e ErfS + J j c dS 

dA A A A 

Charge conservation is ensured by the continuity equation for (p,j c ) or its 
integral equivalent, the Gauss-Ostrogadski law 

(50) J j c rfS = -±-J p dV 

ac c 

valid for arbitrary test volumes C. Of course Faraday's law (Maxwell's 2nd 
integral equation) holds unalteredly 

(51) - J Eds = J M^dS + J K ro HdS 

dA A A 

Perhaps, in the present context the loss currents called forth by the electric 
and magnetic conductivities n e , K m in Maxwell's equations (|4"5)l . (|5T|) attract 
some attention. They simply allow for simulating dissipative effects of any 
origin in the plasma. Beyond it, further simplifying assumptions obviously 
enter the above description, which may be modified or dropped in each case. 
Thus, the model neglects the direct repulsive Coulomb interaction between 
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the particles, which amounts to a low density approximation. Also, the 
fixed mean collision frequency v c may be substituted by a state dependent 
expression. 

In the TLM model the physical state, given by the Maxwell field (E, H), 
the charge density p, and the mean particle velocity v has to be related 
to cell states in a suitable interpretation. In a convex hexahedral parcel 
twines cell (wherein ports on the cell faces interconnect the midpoints of 
the cell edges, cf. fig.l) the fields are connected to cell states in terms of 
port voltages and currents 

(52) f/ M :=E.p^ , := +(-)H • p" 
the sign corresponding with that in 

(53) f=+(-)p"Ap y 

where f denotes the cell face vector which points into the cell. On the cell 
surface t/ p , K as finite integrals approximate line integrals 

(54) UP a J Eds ; I* a + (-) j Hds 

(every port vector in parcel-twines location is naturally identified with an 
integration path on the pertinent face). 

J7™, 7" measure the field strengths in the cell according to (|52|) . yet with 
the nodal images p M,I/ , cf. sect. 2, (HJ[2]). The model equations for Ampere's 
and Faraday's laws are developed in detail in references [T^], [2], [TS]. A 
suitable change of representation, cf.[15] sect. 4, provides decoupled equa- 
tions in terms of transformed voltages and currents (w^' p , ^' p ) M =i....,i2- On 
the first components z = (1x^,1^)^=1 3 the discretized 1st Maxwell equa- 
tion reduces to 

(55) it) = f z n ( t + T - ) + Vl z n { t - T - ) + J c 

the convection current J c being added here as the following perturbation 

(56) J c = I d et(5) J B- 1 |v= \b- x Qv 

Q is the total charge in the cell and v the mean particle velocity. The node 
vector matrix B with adjoint is defined in reference [15] and (modulo a 
factor 4) also in reference |14j . In the present normalization the determinant 
of B equals the cell volume and A = det(B)B~ 1 is the so called area matrix, 
cf.[15] sect. 3. 

Still Q = pV and v have to be incorporated into the TLM model as cell state 
functions. This is achieved by introducing six additional ports in the form 
of the cell face vectors f M which measure the convection currents through 
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the respective cell boundary faces and through pertinent image areas in the 
node. The link states associated with these ports are 



(57) J M 



Jin ~\~ J out 
Jin J out 

to which the following interpretation is construed 

(58) J u := J^l = J a ,vn + Jn.out = PV ■ f M 

Note that voltages are not defined in parallel with these currents. This 
amounts to neglect a pressure in the here adopted low density i approxima- 
tion. In transmission line parlance, the characteristic impedance i i of the 
cell face channels is zero, just as then are the voltages U n = lJ n ,2, to which 
a pressure p — pU is assigned. The generalized TLM scheme, which in the 
setting of sect. 3 dispenses with impedances, still works for 3 = with link 
states defined as in (57). The total charge Q in a cell is related to the cell 
boundary current by charge conservation cf. (|47|) , which yields the updating 
instructions 6 

(59) Q(t + r)=Q(t)+rY / JP(t) 

In the cited canonical representation, [15] sect. 4, the nodal states z% = 
( u fc'*fc) are interpreted as field strengths 

(60) E = B~ x (ufc) fe=li2 ,3 5 H = B' 1 (<JJ) fc=7i8)9 

These substituted in equation (|44[) yield the Lorentz force and updated 
particle velocities in the cells 

(61) v(t + r) =v(t) + TM- 1 F L 

for M(v), Fl(v, E, H), and v at time t. 

Ultimately, equations ([55|) , (|58|) , (f59|) , (|6Tj) are simultaneously solved by 
deflection, branching out the trivial process J n (t) = and the reflection 
map 31 of the free Maxwell field model as outlined in sect. 4. This results in 
the following updating instructions for the coupled TLM process 

u=l 

-It 



(62) 


z out {t 


Hr) 


(63) 




Q 


(64) 


v(t - 




(65) 




Jc 


(66) 
(67) 


T)(t- 

J^L,OUt " 


hr) 
Hr) 


(68) 


Q(t- 





v(t) + tM Fl 



(TToutO^o) 1 {(^0 - <Pl)D(t) + J c } 



6 



= q + r) J? t (t + T) 



/ j J fl,OUt \ 
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with all operations carried out in this order. 

Stability requires the assignments must be contractive in the domain of 
operation ( i.e. up to any allowed maximum excitation). This is ensured 
with bounds for the time step following the guide-lines of sect. 6. 

A. Appendix 

To remain technically in contact with the familiar framework (which uses 
line impedances) the bridge is thrown in this passage from Maxwell's equa- 
tions in the transmission line setting [12], [2], [15] to the presently gen- 
eralized formulation that works with projections instead. In the canonical 
representation of the link state space, [15] sect. 3, with 



(69) z = 




(y = 1/3 the characteristic field admittance) the generalized Ampere's law 
discretized in a parcel-twines cell takes the form of equation ([TT]) with 




(70) Vo 
(0 n , l n denote respectively the n x n zero and unit matrix blocks) and 
(71) 





o 3 \ 




r 






\ o 3 


o 3 / 


\ o 3 o 3 / 



with real selfadjoint (i.e. symmetric) 3x3 matrix operators 

(72) T +( _) = \ det(B)i?- 1 + (-)£) B~ l 

B is the node vector matrix, defined ibid., which depends only on the cell 
geometry. n e and e denote respectively the electric conductivity and per- 
mittivity tensors. With (|69p the projections into the incoming and outgoing 
states related to Ampere's law (separated) are 

(n\ 1 I 13 313 I 

(to) IT i n — - , TT out 



2 



2/13 1 



3 




Proceeding as outlined in sect. 5, i.e. substituting for z n ' p in equation pip 
with operator coefficients (|70p , ([7T]) the total scattering response of a Dirac 
pulse, cf. (|2"6"]l . (|2"T|) . yields the S-matrix blocks K^, . . . , Na pertinent to 
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the generalized Ampere's law. Modulo gauge transformations (|26|) they are 
uniquely determined as 



(74) 

with 
(75) 



Ka = TTout U n in 
Na = n s Vir s 
L A =ir out W(ld-N 



2\3 



u 



M A = n s (Id-Niyw*** 



V 



VT+ 1 - 1 3 


o 3 


3 


-yT+ 1 + 1 3 



yTi 1 - T+ X T_ 



W maps the stub vectors pertinent to Ampere's law back into the link 
channels and acts in the canonical representation on zf © z s as a matrix 
operator 



(76) 



W = ■KfW'Ks = 2 




It is left as an exercise to show that S7 with these data is unitary iff n e — 0. 

For Faraday's law the derivation of operators Kp . . . Np runs completely 
parallel. 



2. Conclusions 



Despite the rich potentiality they offer, the gauge properties of the TLM 
scheme remain still largely unexploited, today, in numerical model design. 
A general method that allows for incorporating a linear or non-linear inter- 
action into a given linear model using these properties has been developed 
in this paper, on the algorithmic level and within a prototype application. 

From a fundamental point of view, future work should yield still deeper 
insight into the precise nature and the range of the observed gauge similarity. 
It is clear from the preceeding that the latter is not to be (mis)understood 
in the restricted sense of Lagrangian field theory: Until now, a Lagrangian 
is not known in the TLM context (though it may be challenging to look 
for such an object). Instead, the structural relationship traced out points 
to some more general common aspects: Much alike gauge field models, the 
TLM algorithm works with unobservable quantities that intermediate inter- 
actions and which together with their internal symmetries can be similarly 
described in fibrations. There is a considerable differential geometric overlap 
with gauge field theory that attracts our attention and further study in this 
direction may significantly promote the TLM method. 

Ultimately, not a few technical spade-work of course remains to be done 
in any particular application that is conceivable. The propagator approach 
marks only a canonical path equipped with reliable instruments. 
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